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We use the non-equilibrium Green function (NEGF) method in the balhstic hmit to provide 
a quantitative description of the conductance of graphene pn junctions - an important building 
block for graphene electronics devices. In this paper, recent experiments on graphene junctions 
are explained by a ballistic transport model, but only if the finite junction transition width, D^, 
is accounted for. In particular, the experimentally observed anamolous increase in the resistance 
asymmetry between nn and np junctions under low source/drain charge density conditions is also 
quantitatively captured by our model. In light of the requirement for sharp junctions in applications 
such as electron focusing, we also examine the pn junction conductance in the regime where Dw is 
small and find that wavefunction mismatch (so-called pseudo-spin) plays a major role in sharp pn 
junctions. 

PACS numbers: 



I. INTRODUCTION 

Graphene is a two-dimensional sheet of carbon atoms 
arranged in a honeycomb lattice with unique electronic 
properties due to its linear energy dispersion, with zero 
bandgap, and a spinor-like two-component wavefunction 
[U 121 [3]. These characteristics give rise to interesting 
transport phenomena such as the absence of backscatter- 
ing |4] , anomalies in the quantum Hall regime |2l [5] , weak 
anti- localization [6l |7] , so-called Klein tunneling [8] , and 
electron focusing analogous to optical effects that occur 
in negative refractive index materials [9 . As such, one 
expects that grahene pn junctions should differ from tra- 
ditional semiconductor pn junctions. The pn junction is 
a basic building block for electronic devices. Developing 
a quantitative understanding of graphene pn junctions is 
an important step on the way to realizing novel devices 
such as graphene lenses [9l [10] and filters [8^ . Our goal 
in this paper is to demonstrate that the non-equilibrium 
Green's function (NEGF) approach [11] quantitatively 
explains recent experiments on graphene pn junctions, 
including the critical role of the junction depletion width 
and the increased in the odd resistance observed under 
low source/drain charge density conditions [T2] . 

Electron transmission across a graphene pn junction 
occurs by interband tunneling. A theoretical treatment 
for an abrupt, graphene pn junction predicts that for a 
symmetric pn junction (i.e. one in which the hole and 
electron densities on each side of the junction are equal) 
the transmission probability is given by cos^O, where 
is the angle between the electron's wave-vector and the 
normal to the junction interface [8 . Realistic pn junc- 
tions will have a transition region of finite width. For 
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a smooth junction transition region of width D^, the 
Wentzel Kramers Brillouin (WKB) approximation can 
be used to show that the transmission probability for 
a symmetric pn junction is ^-^'^fD^sm {e)/2 \^'^ ^ where 
/cj is the Fermi wave- vector. Whether the transition re- 
gion is abrupt or graded, the transmission is perfect when 
^=0 (i.e. commonly refered to as Klein tunneling), but 
the transmission decreases as and Dyj increase. This 
angular selectivity for electron transport across the pn 
junction serves as a filter, allowing states with \0\<(Je 
(where (Jq is the spread of the angular distribution) to 
pass through more effectively. The quantity 2cfq can be 
viewed as the bandwidth of this filter and is what gives 
rise to the larger resistance of a pn junction as compared 
to a uniform graphene sheet. 

Several research groups have recently fabricated 
graphene pnp devices by using electrostatic gates to cre- 
ate p and n regions [121 (HI (HI HI] • The typical setup 
consists of a back-gate and top-gate, which are used to 
control the amount of charge density in the source/drain 
and channel regions respectively. For example, the back 
gates can be set to produce n-type source and drain re- 
gions and the top gate can be biased to change the middle 
(channel) region from n to p type. An asymmetry in the 
device's source to drain resistance as a function of top- 
gate voltage has been experimentally observed [12]. The 
amount of this resistance asymmetry is a measure of the 
intrinsic property of the pn junction, provided that the 
mean free path of carriers is larger than the transition 
length of the junction. One can theoretically compute 
the junction's transition length accounting for non-linear 
screening effects [17 . For recent experiments [12 , typical 
transition lengths for pn junction are less than lOOnm. 
Recent experiments indicate that the carrier's mean free 
path is about lOOnm under low temperature and moder- 
ate carrier density conditions of IQ-^^cm"^ [18]. In addi- 
tion, there is experimental evidence of Fabry- Perot inter- 
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ference effects within the channel region in devices with 
channel lengths less than lOOnm [HI [20], evidently a sig- 
nature of phase coherent transport. Therefore, a ballistic 
transport model is sufficient for the study of the experi- 
mental pn junction devices reported in [12 . 

In this paper, we present a quantitative study of the 
near-equilibrium IV characteristics of graphene pn junc- 
tions. In particular, we present a systematic study of the 
impact of the junction transition width, D^, on the trans- 
port properties of graphene pn junctions. We use the 
non-equilibrium Green function (NEGF) approach with 
a nearest neighbor tight-binding description of graphene. 
This method allows us to accurately simulate the pn junc- 
tion conductance for both abrupt and graded junctions 
under different bias conditions. The value of is not 
known a priori because it depends on charge screening 
and the gate potential as governed by the Poisson equa- 
tion. In this work, we employ the analytical screening 
model derived in [17 . Our numerical result based on the 
assumption of ballistic transport is in reasonable quanti- 
tative agreement with the experiments reported in [12 . 
As has been noted previously, we also found that the 
conductance follows an inverse square root dependence 
on when is large [13], but we find strong de- 
viations from this trend occur when is sufficiently 
small. Understanding this regime of operation has prac- 
tical importance because devices based on the electron 
focusing property of graphene pn junctions, i.e. graphene 
lenses [HITO] and filters [8 , are expected to operate in this 
regime. 



II. GRAPHENE PN JUNCTIONS AND 
RESISTANCE ASYMMETRY 

Before describing the simulation method used in this 
study, we define some terms that will be used in subse- 
quent discussions, present a simple picture of the conduc- 
tance of a graphene npn or pnp junction, and identify the 
issues that will be addressed in the numerical study. Fig- 
ure [TJa,b) are depictions of the energy band diagram for 
the experiment, which shows the location of the intersec- 
tion of the conduction and valence bands (the so-called 
Dirac point) vs. position. A back gate controls the loca- 
tion of the Fermi level, Ef^ in the source/drain regions. 
For Ef>0 (above the Dirac point), there is an increase 
of electrons with respect to Ef=0 , the electron density 
is greater than zero, ni>0, and the material is n-type. 
If the back gate is biased negative so that £^/<0, then 
ni<0 and the source/drain regions are p-type. Similarly, 
the channel (middle) region can be biased by a top gate 
to be either p-type (n2<0) or n-type (77-2 >0). The top 
gate also controls the built-in potential, V^^, of the np/nn 
junction. Thus, appropriate top and bottom gate volt- 
ages can produce npn, nnn, pnp, or ppp structures. Near 
equilibrium conditions are assumed, so that the source 
and drain Fermi energies are the same. In this paper, 
we shall assume that the channel length is greater than 
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FIG. 1: The energy band diagram of graphene npn and nnn devices 
depicted in (a) and (b) respectively, (c) is an intensity plot of the 
pn junction conductance, CTpn, as a function of Ef and Vpn using 
Eq. [l] and (d) plots apn as a function of Vpn at Ef =0.3eV using 
Eq. 



the carrier's phase coherent length, allowing us to ignore 
wave interference effects within the channel. Therefore, 
we can reduce the problem by only studying the trans- 
port across a single pn junction. 

As a first step to understanding the near-equilibrium 
conductance of a pn junction, apn (S/um), as a function 
of Vpn and Ef conditions, we adopt a simple density-of- 
states argument in the Landauer picture. In this simple 
analysis, (jpn can be written in the following form. 



2e' 

7 

Wh 



lin (Mi,M2) 



(1) 



where M1/2 is the number of modes in the source/channel 
respectively and W is the device width. Eq. [l] mim- 
icks the matching of transverse momentum between the 
source and channel in a ballistic manner, such that the 
current will be limited by the region with the smaller 
number of modes. In this simple exercise, we assume a 
zero temperature treatment, so that M1/2 refers to the 
number of modes at their respective Ef. 

Figure Tic) is an intensity plot of apn computed us- 
ing Eq. fTTversus Vpn and Ef. The dark blue lines are 
regions of low current. The horizontal line occurs when 
the source Fermi level is at the Dirac point. Since ni=0, 
so no current fiows. The diagonal black line described by 
Vpn=Ef occurs when the channel region is adjusted to 
place the Dirac point in the channel at the Fermi energy 
of the carrier. Since n2=0 along this diagonal line, no 
current flows. Very similar features are observed experi- 
mentally [12 , 19 . Looking more closely, we can plot apn 
vs. Vpn at a fixed Ef=0.3eV. As shown in Fig. [l](d), the 
conductance vs. Vpn is symmetrical about Vpn=Ef. In 
experiments, an asymmetry about Vpn=Ef is observed 
[12] H^. This asymmetry cannot be captured by sim- 
ple density-of-states argument because its origin is quan- 
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FIG. 2: Total transmission as a function of energy calculated for 
(a) 55 layers (W ^ 13.7nm) metallic armchair ribbons and (b) 32 
layers (W 8nm) semiconducting armchair ribbons transporting 
across a pn junction with a built-in potential of Vpn=0-8eV . In 
both cases, we compared the method outlined in our paper (solid 
line) with that of real-space tight-binding approach (open symbols). 



turn mechanical in nature i.e. quantum tunneling and 
wavefunction mismatch. Our goal in this paper is to ex- 
plain the role that quantum tunneling and wavefunction 
mismatch played in the observed asymmetry and then 
quantitatively explain the magnitude of the asymmetry 
observed in experiments. 



III. SIMULATION METHODS 

In this section, we briefly outline the NEGF formal- 
ism for quantum ballistic transport and describe its ap- 
plication to non-equilibrium transport [HI |2T] through 
graphene pn junctions ^36^. The central quantity in the 
theory is the retarded Green's function, 



G{e, ky) = [(e + i6)I - H{ky) -U - S^(e, ky 



where 5 is an infinitesimal quantity in the channel region 
but is adjusted to provide a non- vanishing density-of- 
states at the Dirac point for the contact regions [22] [23l. 
The Hamiltonian is formulated by treating only nearest- 
neighbor interaction between the Pz orbitals [24l [25] . As- 
suming that the device width is large and homogeneous 
along the direction transverse to current flow, we can 
write the Hamiltonian as 
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where a. Pi and P2 are all 2 x 2 matrices given by. 
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where tc is the nearest neighbor orbital coupling energy 
and ty=tc -\- tce'^^y^^ . The lattice parameter, ao=\^acc 



, where acc=^-^^^ is the c-c bond distance. The quan- 
tum number ky is the quantized transverse momentum, 
to be elaborated upon in next paragraph. The contacts' 
self-energies, Hi(e, ky)^ are obtained by solving 12i=rigiTj 
where gi is the surface green function associated with the 
contact. We should mention that for armchair ribbon an 
analytical closed form solution for gi is possible, since the 
wavefunction is known, both in the tight-binding [26] and 
Dirac formalisms [2T. However, in our numerical treat- 
ment, we had employed a non- negligible 5 in the contact 
regions, which prevent us from using the analytical closed 
form solution for gi. Therefore, gi is computed numeri- 
cally with an iterative scheme proposed in [28]. Finally, 
the current through contact i can then be computed us- 
ing. 
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where A=i{G — G^) is the local density-of-states, 
I]*^=/o(e)r^(e) is the filling function (analogous to the 
in-scattering function for incoherent case), /o(e) is the 
Fermi function of the contacts, and r^=i(I]^ — Sj) is the 
contact broadening factor. In Eq. [sj G^(e) is the elec- 
tron correlation function given by G(I]^^ + E*^)G^. G 
and G^ are computed using the recursive green function 
algorithm through Dyson's equation [29], and making use 
of the fact that the Hamiltonian is tridiagonal in nature. 

By assuming an armchair ribbon configuration and im- 
posing a box-boundary condition and solving the Dirac 
equation, Brey and Fertig showed [2T that the transverse 
momentum is quantized according to 
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for all integer n and W is the width of the device. The 
last term accounts for the momentum of the Dirac points, 
K and K' ^ where the upper /lower sign is used when n is 
even/odd respectively. 

Treating the problem in terms of transverse modes 
greatly reduces the computation burden while still pro- 
viding accurate results when the potential in the trans- 
verse direction is uniform [30l [STJ [32]. This approach 
essentially translates a two-dimensional real space trans- 
port problem into m decoupled one-dimensional real 
space transport problem, where m is the number of rele- 
vant transverse modes. In the limit of large device width, 
VK, we have mocWe. Fig. [2] compares this mode space 
method with two-dimensional real space NEGF calcu- 
lations for a 55-layer metallic graphene ribbon and a 
32-layer semiconducting ribbon, 'n-layer' refers to the 
number of layers of carbon atoms along the width direc- 
tion, where W = na^. In both cases, the two methods 
give nearly identical energy-resolved transmission func- 
tions for transport across a pn junction. The mode space 
approach has a computational burden that scales linearly 
with device width, which makes it possible to study typ- 
ical experimental npn-type structures. 
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FIG. 3: Simulation results from NEGF calculation of an abrupt 
graphene pn junction with device width of 0.5um at temperature 
of 4:K. (a) shows the intensity plot of conductance as function of 
built-in potential Vpn and Fermi energy Ef (see inset of figure |2|b) 
for definitions). Blue/red denotes low/high intensity respectively, 
(b) and (c) plots the conductance curves as function of Vpn (at 
Ef=0.3eV) and Ef (at Vpn=0-QeV) respectively. The dashed lines 
indicates the estimated conductance from a simple density-of-states 
argument (see text) for each case. In each case, the conductance of 
an unbiased graphene (dnn) and a symmetric pn junction (cTpn) are 
indicated, (d) depicts the constant energy contour and its pseudo- 
spin alignment at each side of the junction for an unbiased graphene 
and a symmetric pn junction. 



IV. RESULTS 



As discussed in Sec. II, the experimental setup 
[121 [111 [ISl US] uses a back-gate and top-gate to con- 
trol the source/drain contacts Fermi energy, Ef^ and the 
potential difference across the pn junction, Vpn^ respec- 
tively. The electron density in the source/drain regions is 
rii and in the channel, n2; a negative value can be inter- 
preted as a positive hole density. In this study, we assume 
that the applied source-drain bias is small, in accordance 
with experiments [12]. 

We shall also assume that the channel length is greater 
than the carrier's phase coherent length, allowing us to 
ignore any quantum interference effects within the chan- 
nel i.e. Fabry Perot effects. On this premise, the trans- 
port process across a pn junction on one side of the chan- 
nel would not be influenced by the presence of the pn 
junction on the other side of the channel. Hence, we shall 
focus only on the transport across a single pn junction. 
We begin by examining the pn junction's conductance 
(cr in units of S/um) when the junction is abrupt i.e. 



A. Conductance for an Abrupt PN Junction 

Fig. |3|a), which shows the computed conductance 
as function of Ef and Vpn for an abrupt pn junction, 
exhibits the features typically observed in experiments 
[T2l [T9] . The four distinctive regime of operations (i.e. 
nn, pp, np and pn) are partitioned by the Ef=0 and 
Ef=Vpn lines, which correspond to conditions where ei- 
ther/both the contacts and channel are at zero equilib- 
rium charge density as previously discussed in Se. II us- 
ing simple density-of-states argument. Examining more 
closely. Fig. |3|b,c) plot the conductance as a function of 
Vpn {Ef=0.3eV) and Ef {Vpn=0.6eV) respectively. The 
red solid line is NEGF calculation while the dashed black 
line is calculated using Eq. [l] In Fig. [sj^b), one can 
clearly see the conductance asymmetry with respect to 
the Vpn=Ef point. This conductance asymmetry is an 
experimentally observed phenomena which cannot 
be captured by a simple density-of-states argument. 

Conductance asymmetry refers to the difference in con- 
ductance between a np junction and its nn counterpart. 
By 'counterpart', we mean that n2 for np junction is 
—712 for nn junction, while ni is the same for both junc- 
tions. In Fig. [sjb), (Jnn is the conductance of an unbiased 
graphene (i.e. Vpn = 0), while apn is the conductance of 
a symmetric np junction (i.e. Vpn = 2£^/). Clearly ann> 
(jpn^ a feature that is observed experimentally and cap- 
tured by our NEGF calculation. The asymmetry is due 
to wave- function mismatch at the junction interface. 

Electrons in graphene have a two-component wave 
function, which is often refered to as pseudo-spin in anal- 
ogy to the two-component wave function for spins. In the 
vicinity of the Dirac point, the two-dimensional Hamilto- 
nian can be written in a form of Dirac equation O [4] i.e. 
H = VfG • where a and P are the Pauli spin matri- 
ces and momentum operator respectively. By convention, 
the definition of pseudo-spin is such that its direction is 
parallel to the group velocity, since the group velocity 
operator is defined as vq = V pH = VfB. Fig. |3|d) 
provides an illustration of the constant energy contour 
for an unbiased nn graphene and symmetric pn junction 
case. The arrows simultaneously represents the group 
velocity and pseudospin, which points inwards/outwards 
for the valence/conduction band respectively. For the 
nn case, the velocity vectors are similarly aligned for 
each side of the junction. For the pn case, vq changes 
sign across the junction. Analogous to spin, the wave- 
function for the n and p regions can be expressed as 
\^n) = (l,e^^)/V2 and = (1, e^(^+^))/y2 respec- 

tively (where = tan~^{ky/kf)). The tranmision proba- 
bility across the junction for a particular mode can then 
be written simply as For ky=0, the wave- 

function is perfectly matched i.e. 1^^) = so trans- 

mission is unity (i.e. Klein tunneling [8j). Through the- 
oretical analysis, we derive apn = 2/3ann^ where the fac- 
tor of 2/3 is due to T^ky{^ — ky/k'j). Similarly, it can also 
be shown that the conductance when Vpn oo would 
approach the asymptotic value of ~(l/2 + 7r/8)(Tnn- Rig- 



5 




0- 7 1 I 1.5 




Transport Direction (nm) Junction Length (nm) Fermi Energy (eV) 



FIG. 4: Analysis of the impact of junction width on conductance modulation of a symmetric pn junction, (a) plots the 
energy-resolved local density- of- states (i.e. for a W=0.5um pn junction device with a junction length of Du,=^Onm and a 

built-in potential Vpn=0-5eV . This is generated for the transverse mode with a transverse energy of ey^0.075eV , which yields an apparent 
'bandgap' of ~0.15ey as seen in the plot, e = is set at the Dirac point of the source (b) is the intensity plot of the transmission as a 
function of incident angle and Dw The device is a symmetric pn junction with Vpn='2E f , where E f=0.3eV . An illustration of the setup 
is shown in the inset, (c) plots the fraction of conductance contribution due to wavefunction mismatch, i.e. cr^rif ^ /cpn , for a symmetric 
pn junction (see Eq. [7|. 



orous NEGF calculations faithfully reproduces these fea- 
tures as shown in Fig. [sjb). 

In summary, the conductance of an abrupt pn junc- 
tion can be understood as being controlled by the region 
in which the number of conducting channels is smallest. 
Wavefunction mismatch reduces the current by a factor 
of 2/3 in a symmetric np junction compared to an un- 
biased graphene. These features are captured by NEGF 
simulation, but realistic junctions have a finite transition 
width, D^, which, as we emphasize in the following sec- 
tion, plays an important fole. 



B. Effect of Junction Width on Symmetric PN 
Junction Conductance 

For realistic np junctions, the transition length across 
which the charge density changes monotonically from n- 
type to p-type is finite. The width of the junction transi- 
tion region has a strong influence on the conductance of 
the junction. To understand these effects, it is instruc- 
tive to consider each of the transverse modes individually. 
According to Eq. [6j for a wide graphene sheet, the trans- 
verse momentum with respect to the Dirac point is given 
hy ky ^ nir /W where n is an integer. We can view each 
mode as a ray with an angle of incidence on the junction 
= tan~^[ky{k'j - ky)~^-^]. Fig. [ij^a) plots the energy- 
resolved local density-of-states (GS*^G^) for a pn junc- 
tion with a transition length of Dyj=50nm and a built-in 
potential Vpn=0.5eV. Only one transverse mode with a 
transverse energy of hvfky^O. 075eV is considered. The 
dark blue region correspond to a low density of states. 
Because of the quantization of transverse wavevectors, 
an apparent bandgap of O.lbeV is observed. There- 
fore, electrons at Fermi energy of Ef<Vpn will have to 
quantum mechanically tunnel through this /c^y-dependent 



bandgap when moving from one side of the junction to 
the other. This problem is analogous to the classic band- 
to-band tunneling problem in direct band-gap semicon- 
ductors [33 . 

Based on this physical picture, Cheianov and Falko 
[13] worked out the WKB tunneling probability for a 
given transverse mode to be g-^^/^^^s^n (6')/2^ This tun- 
neling expression is derived by assuming a symmetric pn 
junction. In similar spirit to classic band-to-band tun- 
neling treatment, the electric field across the pn junc- 
tion is assumed to be linear, given by 2Ef/D^ in this 
work. Realistically, the potential energy landscape at 
the beginning/end of the pn junction would exhibits a 
more quadratic profile. However, only the details of 
the potential energy landscape within the tunneling re- 
gion contributes to the WKB tunneling probability. In 
this linear electric field approximation, the tunneling dis- 
tance is simply given by, Dt = 2hvfkyD^/Vpn' Fig. 
|4jb) shows the NEGF-computed electron transmission 
as a function of incident angle, ^, and transition length, 
D^, for a symmetric pn junction with Ef=0.3eV. As 
expected, increasing results in a decreased angu- 
lar bandwidth (2<J6>) of the allowable transverse modes 
that can be transmitted across the pn junction, which 
subsequently leads to a decreased pn junction conduc- 
tance. Based on the WKB tunneling forumla, it can be 
shown that Gpn^^JkjjD^. This leads to the question 
of whether the junction rectification metric can be im- 
proved by using a larger D^. Unfortunately, (7nnl^T>n 
which depends on Dyj in an inverse square root manner, 
yields only a factor of 10 with Dyj ^ lOOnm at a typical 
Ef of 0.3eV. 

Device concepts based on the electron focusing prop- 
erty of pn junction, i.e. graphene lenses [9 , 10 and filters 
[8 operate best in a symmetric pn junction. These de- 
vices operate in a regime where (Jq has to be as large as 
possible so as to reconstruct back a point source image 
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on the other side of the pn junction. This imphes that 
Dyj has to be sufficiently smah to enhance tunnehng at 
high incident angle. By accounting for both the wave- 
function mismatch and tunneling factor in the following 
manner, T{0) = cos2(l9)e-^^/^-^^^'^^)/^ we found that 
we are able reproduce the result of NEGF for arbitrary 
(not shown). By integrating T{0) over all transverse 
modes, we can arrive at a more general result for the con- 
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ductance of a symmetric pn junction (in units of 2^), 



-^Dqj^ k j 



V2 



:erf{^jDu,kf) 
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where 7 = 7r/2. The first term in Eq. [t] is due to 
the tunneling factor. The last two terms are corrections 
due to wavefunction mismatch. Eq. [7| can be written 
as apn = (Jp^ + ^pn^' Fig. [ij^c) is an intensity plot 
of cFpi^ /cTpn as a function of Ef and Du,. The blue 
region represents ^pn^ while the red regions 

indicate <C o-pn^- Evidently, conductance modula- 
tion is predominantly due to wavefunction mismatch only 
when ^D^kf < 1, which suggests that for the tunneling 
component not to limit electron focusing applications, a 
< 5nm is required. 
In summary, the conductance asymmetry of a graphene 
pn junction is due to two quantum mechanical processes, 
wavefunction mismatch and the need to tunnel through 
an apparent bandgap induced by the quantization of 
transverse momentum. Increasing results in a de- 
creased angular bandwidth {2ae) of the allowable trans- 
verse modes that can be transmitted across the pn junc- 
tion. This leads to a decreased pn junction conductance 
and would eventually result in a larger magnitude of odd 
resistance Rodd- 



C. Odd Resistance and Comparison With 
Experiments 

Finally, we shall examine the odd resistance, Rodd^ 
of pn junction devices and compare our NEGF result 
with the experimental data reported in [12]. Typically, 
the resistance asymmetry is characterized by analyzing 
the resistance of the device as a function of Vpn at a 
given Ef. A quantity known as odd resistance, Rodd^ 
can be obtained by taking the difference between the re- 
sistance of the npn and its nnn counterpart device i.e. 



R. 



odd 



■npn 



Rr) 



where the channel hole density 



for npn is equal to the electron density for nnn. Fig. 
[TJa,b) depicts the energy band of a typical npn device 
and its nnn counterpart. Rodd for a long channel de- 
vice depends only on the odd resistance contribution from 
individual pn junctions. Essentially, Rodd simply reduces 
to Rodd = Wl^np — ^/(^nn] obtained by a simple sum of 
the resistance of the two adjacent pn junctions. As a 
ffist step, we shall investigate the contributions of wave- 
function mismatch and quantum mechanical tunneling 
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FIG. 5: Theoretical resistance per unit width of a single pn junc- 
tion as a function of Vpn at Ef=Q.lSA^eV under different D uj con- 
ditions (b) Comparison of NEGF and experimental odd resistance 
I^odd for different Dw plotting with respect to 722 — ni. Experi- 
mental device has W=1.7um with contacts Fermi energy also at 



=0.1845ey. 
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FIG. 6: Comparison of the NEGF simulated Rodd with that 
of experimental results in 12 at different Fermi energies i.e. 
£;/=0.1917, 0.165, 0.1381 and' 0.1044ey. The device width is 
1.7um and operating temperature is 4K. In the simulations, the 
contacts are assumed to have an energy broadening of ImeV. The 
calculations assumed q;=0.78 due to a of 4.5 12 and an oxide 
thickness of Tocc~80nm as a fitting parameter. 



processes to the magnitude of Rodd- Fig. |5|a) shows the 
theoretical pn junction resistance as a function of Vpn at 
Ef=0.184:5eV under different D^, conditions. The odd 
resistance contribution due to wavefunction mismatch 
alone (i.e. D^=0) does not adequately account for the 
Rodd observed in experiments as evident in Fig. [5|b). In 
fact, it only accounts for 10% of the Rodd- Accounting 
for finite is essential to match the experimental data. 

In our NEGF simulation, is the parameter that we 
need to determine prior to our NEGF calculations. The 
determination of is an electrostatics problem which 
is sensitive to the specific device geometry. In this work, 
we used the model presented in [17 , which expresses 
in terms of ni and 712 as follows. 



: 0.196 X 



hvfa 



1/3 



n2 
ni 



722 



(8) 



where a=e^/ (tvhvf). tz and Tox are the effective dielectric 
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constant and thickness of the oxide between the top gate 
and graphene device respectively [37 . 

Fig. |6] shows the computed Rodd at different Ef con- 
ditions as a function of ni — 712 and its comparison with 
experimental data. The NEGF result achieves quanti- 
tative agreement with the experimentally observed odd 
resistance. In particular, the increase in odd resistance 
with decreasing a puzzling feature in the experiments 
[12], is captured by the simulations. This occurs because 
of the increasing with decreasing Ef i.e. smaller 
Til which results in an increase of pn junction screening 
length. We note conductance oscillations at the small- 
est Ef. These oscillations are likely due to interference 
effects within the device channel, which are more pro- 
nounce at Ef =0.1eV due to the large which leads 
to an effectively shorter channel length. The 'spikes' ob- 
served in the NEGF simulations when Rodd crosses zero 
are due to the zero density-of- states at Dirac point. Such 
spikes are not observed in the experiments, probably due 
to the presence of spatial fluctuations (electron-hole pud- 
dles) when Ef approaches the Dirac point [34^. By con- 
struction, our NEGF model does not account for these 
electron/hole puddles. 

V. CONCLUSIONS 

In this paper we presented a numerical study of elec- 
tron transport in graphene pn junctions. We first pre- 



sented a very simple minimum density-of-states (or con- 
ducting channels) model to account for the overall shape 
of the conductance vs. source carrier density and junc- 
tion potential. Such a simple model does not capture the 
resistance asymmetry observed experimentally. We then 
use NEGF simulation to explore in detail the role of wave- 
function mismatch (also called pseudo-spin) and quan- 
tum mechanical tunneling through the junction transi- 
tion region. In particular, we examined deviations from 
the inverse square root dependence of apn on D^j due to 
wave-function mismatch at small D^. Finally, we com- 
pared the simulations to a recent experiment and showed 
that the numerical model is in reasonable agreement with 
experiments and explain the increase in odd resistance 
with decreasing carrier density in the source. The novel 
features of graphene's electronic structure lead to inter- 
esting possibilities for new devices, and this study shows 
that NEGF simulation should provide a useful tool to 
explore and assess device concepts. 
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limit where it is equivalent to the Landauer approach, as 
explained in chapter 8 of [11] . 

We shall also briefly discuss the validity of the analytical 
non-linear screening model of Eq.js] As discussed in [17] , 
the accuracy of the analytical model is dependent on the 
physical parameter a. The model would yield a possible 
deviation from numerical results of ~ 25% when a ~ 0.9 
at the charge neutrality point within the transition re- 
gion. We emphasize that the error is of an oscillatory 
nature with the oscillation amplitude deterioting when it 



moves away from the charge neutrality point. This would 
serve to alleviate the average error. While on the other 
limit, a ~ 0.1 yields an excellent agreement with the nu- 
merical result. In the set of experiments we studied in this 
manuscript, a ~ 0.78 due to a At ~ 4.5. Therefore, the av- 
erage error is intermediate of these two limits. However, 
when a is large (approaching 1), one would need to also 
account for electron exchange and correlation effects (see 
e.g. 35 ), making the problem numerically non-tractable. 



